
#include "lagrange_1d.hh"

#include <iostream>
#include <iomanip>
using namespace std;

template<unsigned int ORDER, typename T>
Point<SPATIAL_DIMENSION, T> LagrangeSegmentFiniteElement<ORDER, T>::vertices[LagrangeSegmentFiniteElement<ORDER, T>::Traits::N_Nodes];

template<unsigned int ORDER, typename T>
T LagrangeSegmentFiniteElement<ORDER, T>::weights[LagrangeSegmentFiniteElement<ORDER, T>::Traits::N_Nodes];

template<unsigned int ORDER, typename T>
bool LagrangeSegmentFiniteElement<ORDER, T>::compute_vertices_called = false;

template <unsigned int ORDER, typename T>
int LagrangeSegmentFiniteElement<ORDER, T>::test ( )
{
  int i, j, n;
  int test_num = 5;

  std::vector<Point<SPATIAL_DIMENSION, T> > nodes ;
  this->reference_element_coords(nodes);

  cout << "\n\nBASIS_LAGRANGE_SEGMENT_TEST: Order=" << ORDER << "\n";
  cout << "  Verify basis functions for 1-D elements.\n";

  n = nodes.size();
  cout << "\n   Number of reference element nodes:= " << n ;
  for (i=0; i < n; ++i)
    cout << "\n        nodes[" << i << "] = [" << nodes[i][0]  << ", " << nodes[i][1]  << ", " << nodes[i][2] << "]" ;
  cout << endl ;

  cout << "\n  The basis function values at basis nodes\n";
  cout << "  should form the identity matrix.\n\n";

  for ( j = 0; j < n; j++ )
  {
    cout << "  ";
    for ( i = 0; i < n; i++ )
    {
      T phi = this->evaluate_basis ( i, nodes[j] );
      cout << setw(12) << phi;
    }
    cout << "\n";
  }
  
  cout << "\n  The X derivatives should sum to 0.\n\n";
  cout << "  \t\tdPhi_i/dX_j ... \t\tsum(dPhi/dX):\n\n";

  // evaluate_gradient (  unsigned int ib, const Point<SPATIAL_DIMENSION>& px, Point<SPATIAL_DIMENSION>& grad)
  Point<SPATIAL_DIMENSION, T> grad ;
  for ( j = 0; j < n; j++ )
  {
    T sum_x = 0.0;
    cout << "  ";
    for ( i = 0; i < n; i++ )
    {
      this->evaluate_gradient(i, nodes[j], grad);
      cout << setw(12) << grad[0];
      sum_x += grad[0];
    }
    cout << setw(15) << sum_x << "\n";
  }
  
  cout << "\n  The X double-derivatives should sum to 0.\n\n";
  cout << "  \t\td^2Phi_i/dX^2_j ... \t\tsum(d^2Phi/dX^2):\n\n";

  // evaluate_gradient (  unsigned int ib, const Point<SPATIAL_DIMENSION>& px, Point<SPATIAL_DIMENSION>& grad)
  PointList<SPATIAL_DIMENSION, T> hess(SPATIAL_DIMENSION) ;
  for ( j = 0; j < n; j++ )
  {
    T sum_x = 0.0;
    cout << "  ";
    for ( i = 0; i < n; i++ )
    {
      this->evaluate_hessian(i, nodes[j], hess);
      cout << setw(12) << hess.x(0);
      sum_x += hess.x(0);
    }
    cout << setw(15) << sum_x << "\n";
  }
  
  PointList<SPATIAL_DIMENSION, T> ips(test_num);

  for ( j = 0; j < test_num; j++ )
    ips.x(j) = generate_random();

  cout << "\n   Number of random evaluation nodes:= " << test_num ;
  for (i=0; i < test_num; ++i)
    cout << "\n        nodes[" << i << "] = [" << ips.x(i)  << ", " << ips.y(i)  << ", " << ips.z(i) << "]" ;
  cout << endl << endl ;

  cout << "\n==   The basis function values at ANY point P should sum to 1   ==\n";
  cout << "   ---------------------P---------------------------     PHI_SUM\n";

  // evaluate_basis_all ( const PointList<SPATIAL_DIMENSION>& ips, std::vector<double>& basis_vals )
  std::vector<T> basis_vals;
  this->evaluate_basis_all/*_unrolled*/ ( ips, basis_vals );

  for ( j = 0; j < test_num; j++ )
  {
    T phi_sum = 0.0;
    for (i=0; i < n; ++i)
    {
      cout << "  " << setw(12) << basis_vals[i+j*n] ;
      phi_sum += basis_vals[i+j*n] ;
    }
    cout << "  " << setw(12) << phi_sum << "\n";
  }

  basis_vals.clear();
  
  cout << "\n==   The basis derivative values at ANY point P should sum to 0   ==\n";
  cout << "   ---------------------P-----------------------     PHI_DERIVATIVE_SUM\n";

  std::vector<Point<SPATIAL_DIMENSION, T> > grad_vals;
  this->evaluate_gradient_all(ips, grad_vals);

  for ( j = 0; j < test_num; j++ )
  {
    T dphi_sum = 0.0 ;
    for (i=0; i < n; ++i)
    {
      cout << "  " << setw(12) << grad_vals[j*n+i][0] ;
      dphi_sum += grad_vals[j*n+i][0] ;
    }
    cout << "  " << setw(12) << dphi_sum << "\n";
  }

  grad_vals.clear();

  cout << "\n==   The basis second-derivative values at ANY point P should sum to 0   ==\n";
  cout << "   ---------------------P-----------------------     DPHI_DERIVATIVE_SUM\n";

  std::vector<PointList<SPATIAL_DIMENSION, T> > hessian_vals;
  this->evaluate_hessian_all(ips, hessian_vals);

  for ( j = 0; j < test_num; j++ )
  {
    T ddphi_sum = 0.0 ;
    for (i=0; i < n; ++i)
    {
      cout << "  " << setw(12) << hessian_vals[j*n+i].x(0) ;
      ddphi_sum += hessian_vals[j*n+i].x(0) ;
    }
    cout << "  " << setw(12) << ddphi_sum << "\n";
  }

  hessian_vals.clear();

  return 0;
}


// instantiations
template class LagrangeSegmentFiniteElement<1, double> ;
template class LagrangeSegmentFiniteElement<2, double> ;
template class LagrangeSegmentFiniteElement<3, double> ;
template class LagrangeSegmentFiniteElement<4, double> ;
template class LagrangeSegmentFiniteElement<5, double> ;
template class LagrangeSegmentFiniteElement<6, double> ;
template class LagrangeSegmentFiniteElement<7, double> ;
template class LagrangeSegmentFiniteElement<8, double> ;
template class LagrangeSegmentFiniteElement<9, double> ;
template class LagrangeSegmentFiniteElement<10, double> ;


